1 Part 1

1.1 Libraries


library(tidyverse)
library(ggspatial)
library(maptools)library(dplyr)
library(ggthemes)
library(spatstat)
library(classInt)
library(readxl)
library(knitr)
library(sf)

1.2 Loading and tidying data

setwd("C:/Users/User/Documents/OneDrive/Desktop/R/spatial")
ge_admin_url <- "https://data.humdata.org/dataset/3ee95199-2dfe-40fc-b9bd-b44cc8c91024/resource/ac2b2747-18a2-43c9-abbb-4ab2c5c539a3/download/geo_adm_geostat_20191018_shp.zip"
download.file(ge_admin_url, destfile="ge_admin.zip")
unzip("ge_admin.zip")
file.remove("ge_admin.zip")
ge_adm_0 <- st_read("geo_admbnda_adm0_geostat_20191018.shp")
## Reading layer `geo_admbnda_adm0_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm0_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 13 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 40.00834 ymin: 41.05455 xmax: 46.72672 ymax: 43.5865
## geographic CRS: WGS 84
ge_adm_1 <- st_read("geo_admbnda_adm1_geostat_20191018.shp")
## Reading layer `geo_admbnda_adm1_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm1_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 13 features and 16 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 40.00834 ymin: 41.05455 xmax: 46.72672 ymax: 43.5865
## geographic CRS: WGS 84
ge_adm_2 <- st_read("geo_admbnda_adm2_geostat_20191018.shp")
## Reading layer `geo_admbnda_adm2_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm2_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 70 features and 19 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 41.54717 ymin: 41.05455 xmax: 46.72672 ymax: 43.25688
## geographic CRS: WGS 84
ge_pop<-read_xlsx("geo_admpop.xlsx",sheet=3)
## New names:
## * `` -> ...8
ge_pop_2 <- ge_pop%>%
  select(ADM2_PCODE, POP_Total,F_Total,M_Total)

ge_adm_2_ii <- ge_adm_2 %>%
  select(ADM2_PCODE, ADM2_EN)

ge_adm_2_pop <- left_join(ge_adm_2_ii, ge_pop_2, by="ADM2_PCODE")

1.3 Adding data for Tbilisi

ge_adm_1_TBI <- ge_adm_1 %>%
  select(ADM1_PCODE, ADM1_EN) %>%
  rename(ADM_PCODE = ADM1_PCODE,
         ADM_EN = ADM1_EN) %>%
  filter(ADM_EN == "Tbilisi")
ge_adm_1_TBI <- ge_adm_1_TBI %>%
  mutate(F_Total=634743,M_Total=536357,POP_Total=sum(634743,536357))

ge_adm_2_pop <- ge_adm_2_pop %>%
  rename(ADM_PCODE = ADM2_PCODE,
         ADM_EN = ADM2_EN)
ge_adm_2_pop <- rbind(ge_adm_2_pop, ge_adm_1_TBI)

1.4 Calculating gender share

ge_adm_2_pop<-ge_adm_2_pop %>% mutate(Share= F_Total / POP_Total)

1.5 Plotting

ggplot()+
  theme_minimal(base_size=8.5)+
  theme(axis.title = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 13,face="bold",hjust=0.5))+
  geom_sf(data = ge_adm_0, fill = "grey", colour = "grey60", size = 0.1)+
  geom_sf(data = ge_adm_2_pop, aes(fill = Share*100),colour = "grey60", size= 0.1)+
scale_fill_gradientn(colors=c("#004648","#E6E2BC","red")) +
  labs(title = "Gender Balance in Georgia",
       fill  = "Share of female \npopulation,%",
       caption = "data: The Humanitarian Data Exchange \nvisual: D.Kajaia")+
  annotation_scale(location = "bl", style = "ticks") +
  annotation_north_arrow(location = "tr", width = unit(0.5, "cm"))

ggsave("Gender_balance_in_Georgia.png", dpi = 600, height = 4, width = 7, units = "in")

1.6 Conclusion

ge_adm_2_pop %>% summarize(sd=sd(Share),mean=mean(Share))
## Simple feature collection with 1 feature and 2 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 41.54717 ymin: 41.05455 xmax: 46.72672 ymax: 43.25688
## geographic CRS: WGS 84
##           sd      mean                       geometry
## 1 0.01318941 0.5123185 POLYGON ((46.39876 41.47027...

From the map we can see that the sex ration is not equal - it is more female-biased. Main reason for that can be that females live longer than males.

As we can see there is some difference in gender share between districts, but Standard deviation of gender share is very small - only 1 %, therefore mostly these differences are uxplainable. It seems that female share is higher in some cities with higher population (for example in Tbilisi). This may be the case due to selective abortion. Selective abortion should be lower in big cities because people in big cities should be more educated.

2 Part 2

2.1 Libraries


library(tidyverse)
library(ggspatial)
library(readxl)
library(knitr)
library(curl)
library(sf)

2.2 Importing data

setwd("C:/Users/User/Documents/OneDrive/Desktop/R/spatial")
url <- "http://aasa.ut.ee/Rspatial/data/FarmedAnimalsByLocation_31102018.xlsx"
destfile <- "FarmedAnimalsByLocation_31102018.xlsx"
curl_download(url, destfile)
agrAnimal <- read_excel(destfile)
download.file("https://geoportaal.maaamet.ee/docs/haldus_asustus/omavalitsus_shp.zip", destfile="omavalitsus_shp.zip")
unzip("omavalitsus_shp.zip")
file.remove("omavalitsus_shp.zip")
list.files(pattern = ".shp")
##  [1] "geo_admbnda_adm0_geostat_20191018.shp"
##  [2] "geo_admbnda_adm0_geostat_20191018.shp.xml"
##  [3] "geo_admbnda_adm1_geostat_20191018.shp"
##  [4] "geo_admbnda_adm1_geostat_20191018.shp.xml"
##  [5] "geo_admbnda_adm2_geostat_20191018.shp"
##  [6] "geo_admbnda_adm2_geostat_20191018.shp.xml"
##  [7] "geo_admbndl_admALL_geostat_itos_20191018.shp"
##  [8] "geo_admbndl_admALL_geostat_itos_20191018.shp.xml"
##  [9] "geo_admbndp_admALL_geostat_itos_20191018.shp"
## [10] "geo_admbndp_admALL_geostat_itos_20191018.shp.xml"
## [11] "gis_osm_buildings_a_free_1.shp"
## [12] "gis_osm_places_free_1.shp"
## [13] "gis_osm_roads_free_1.shp"
## [14] "gis_osm_waterways_free_1.shp"
## [15] "gps_us.shp"
## [16] "omavalitsus_20200504.shp"
municip <-  st_read("omavalitsus_20200504.shp") 
## Reading layer `omavalitsus_20200504' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\omavalitsus_20200504.shp' using driver `ESRI Shapefile'
## Simple feature collection with 79 features and 5 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 369032.1 ymin: 6377141 xmax: 739152.8 ymax: 6634019
## projected CRS:  Estonian Coordinate System of 1997

2.3 Cleaning data

agrAnimal$county <- NULL
agrAnimal$`admin unit` <- NULL
agrAnimal <- agrAnimal %>% rename(action_place = `action place`,
                                  estonian_holstein_cattle = `estonian holstein cattle`,
                                  estonian_red_cattle = `estonian red cattle`,
                                  estonian_native_cattle = `estonian native cattle`,
                                  beef_cattle = `beef cattle`,
                                  y = `X koordinaat`,
                                  x = `Y koordinaat`)
agrAnimal <- agrAnimal %>% mutate(x = as.numeric(x), y = as.numeric(y))
agrAnimal_2 <- gather(agrAnimal, "key", "value", 2:8)
agrAnimal_2 <- agrAnimal_2 %>% filter(value > 0)

agrAnimal_2 <- agrAnimal_2 %>%
  mutate(municip_key = str_to_lower(municipality))

municip <- municip %>%
  mutate(municip_key = str_to_lower(ONIMI))

2.4 Spatial join

agrAnimal_2_sf <- st_as_sf(agrAnimal_2, coords = c("x", "y"), crs = 3301)
agrAnimal_2_sf_municip <- st_join(st_transform(agrAnimal_2_sf, 3301), st_transform(municip, 3301), join = st_intersects)

st_geometry(agrAnimal_2_sf_municip) <- NULL
agrAnimal_2_sf_municip_aggr <- agrAnimal_2_sf_municip %>%
  group_by(OKOOD, key) %>%
  summarise(sum = sum(value)) %>%
  ungroup()
agrAnimal_2_sf_municip_aggr <- agrAnimal_2_sf_municip_aggr %>%
  spread(key, sum)

agrAnimal_2_sf_municip_aggr <- left_join(municip, agrAnimal_2_sf_municip_aggr, by="OKOOD")

2.5 Pigs distribution in municiaplaities

ggplot()+
  theme_grey(base_size=8.5)+
  theme(plot.title = element_text(size = 13,face="bold"))+
  geom_sf(data = agrAnimal_2_sf_municip_aggr, aes(fill=pigs), size=0.25, colour = "grey70")+
  scale_fill_gradientn(colours = c("#B0B579","#5E623F"), na.value = "#B6D3C1")+
  labs(fill = "N",
       title = "Pigs in Estonian Municipalities",
       subtitle = "Agricultural Registers and Information Board",
       caption = "Author: D.kajaia")

2.6 Density

agrAnimal_2_sf_municip_aggr$area <- st_area(agrAnimal_2_sf_municip_aggr)
agrAnimal_2_sf_municip_aggr <- agrAnimal_2_sf_municip_aggr %>%
  mutate(area = as.numeric(area) / 1000000) 
ggplot()+
  theme_grey(base_size=8.5)+
  theme(plot.title = element_text(size = 13,face="bold"))+
  geom_sf(data = agrAnimal_2_sf_municip_aggr, aes(fill = pigs/ area), size=0.25, colour = "grey70")+
  scale_fill_gradientn(colours = c("#B0B579","#5E623F"), na.value = "#B6D3C1")+
  labs(fill = "N per km2",
       title = "Pigs Density in Estonian Municipalities",
       subtitle = "Agricultural Registers and Information Board",
       caption = "Author: D. Kajaia")+
   annotation_scale(location = "bl", style = "ticks") +
   annotation_north_arrow(location = "tr", width = unit(0.5, "cm"))

2.7 Finding municipality

agrAnimal_2_sf_municip_aggr%>%
  filter(pigs==max(pigs,na.rm=TRUE))%>%
  select(ONIMI,OKOOD,pigs)
## Simple feature collection with 1 feature and 3 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 576487.5 ymin: 6437907 xmax: 628026 ymax: 6497708
## projected CRS:  Estonian Coordinate System of 1997
##           ONIMI OKOOD  pigs                       geometry
## 1 Viljandi vald  0899 78237 MULTIPOLYGON (((621049.2 64...
agrAnimal_2_sf_municip_aggr%>%
  mutate(density=pigs/area)%>%
  filter(density==max(pigs/area,na.rm=TRUE))%>%
  select(ONIMI,OKOOD,density)
## Simple feature collection with 1 feature and 3 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 576487.5 ymin: 6437907 xmax: 628026 ymax: 6497708
## projected CRS:  Estonian Coordinate System of 1997
##           ONIMI OKOOD  density                       geometry
## 1 Viljandi vald  0899 57.03777 MULTIPOLYGON (((621049.2 64...

Both – the number and the density of pigs are highest in Viljandi.

3 Part 3

3.1 Libraries


library(rnaturalearth)
library(lubridate)
library(ggspatial)
library(tidyverse)
library(classInt)
library(stringr)
library(readxl)
library(knitr)
library(units)
library(sf)

3.2 Importing data

setwd("C:/Users/User/Documents/OneDrive/Desktop/R/spatial")
dt<-read.csv("https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv")

countries50 <- ne_download(scale = 50, type = 'countries', category = 'cultural', returnclass = "sf")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\User\AppData\Local\Temp\Rtmp0YCM6x", layer: "ne_50m_admin_0_countries"
## with 241 features
## It has 94 fields
## Integer64 fields read as strings:  POP_EST NE_ID

3.3 Cleaning data

recovered <- gather(dt, "key", "value", 5:ncol(dt))
recovered$key<-str_replace_all(recovered$key,'X',"")
recovered<-recovered %>% mutate(place=paste0(Country.Region, ", ", Province.State),date=mdy(key))
recovered$place<-gsub(", $", "",recovered$place)


recovered_sf <- st_as_sf(recovered, coords = c("Long", "Lat"), crs = 4326)
recovered_sf_latest <- recovered_sf %>%
  filter(date == max(date)-days(1) | date == max(date)- days(8))
recovered_crd_latest <- st_coordinates(recovered_sf_latest) %>%
  as_tibble()

recovered_sf_latest_2 <- bind_cols(recovered_sf_latest, recovered_crd_latest)

recovered_sf_latest_2 <- recovered_sf_latest_2 %>%
  filter(X != 0) %>%
  filter(Y != 0)
recovered_sf_latest_3 <- recovered_sf_latest_2 %>%
  select(place, X, Y, value, date) %>%
  st_drop_geometry() %>%
  spread(date, value)

dates <- colnames(recovered_sf_latest_3)[4:5]
colnames(recovered_sf_latest_3)[4:5] <- c("a", "b")
recovered_sf_latest_3$a[is.na(recovered_sf_latest_3$a)] <- 0
recovered_sf_latest_3$b[is.na(recovered_sf_latest_3$b)] <- 0

recovered_sf_latest_3 <- recovered_sf_latest_3 %>%
  mutate(change = b - a)

3.4 Plotting

ggplot()+
  theme_grey(base_size=8.5)+
  theme(axis.title = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 13,face="bold"))+
  geom_sf(data = countries50, colour = "grey", fill= "white", size=0.2)+
  geom_point(data = recovered_sf_latest_3, aes(x = X,
                                               y = Y,
                                               size = change,
                                              colour = change,
                                           alpha= change))+
  scale_alpha(range = c(0.1, 1))+
  scale_colour_gradientn(colours= c("red1", "red4"))+
  scale_size_continuous(range = c(1, 5))+
  labs(size= "Recovered cases", alpha = "Recovered cases", colour = "Recovered cases",
       title = "Recovered from COVID-19 during the last week")+
  guides(alpha = F,
         size= F)

recovered_sf_latest_3_change_top10 <- recovered_sf_latest_3 %>%
  arrange(-change) %>%
  head(n = 10)

ggplot()+
  theme_grey(base_size=8.5)+
  theme(axis.title = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 13,face="bold"),
        plot.subtitle = element_text(size=10))+
  geom_sf(data = countries50, fill = "white", colour = "grey", size=0.2)+
  geom_point(data = recovered_sf_latest_3_change_top10, aes(x = X,
                                               y = Y,
                                               size = change,
                                              alpha = change), colour = "blue")+
  scale_size_continuous( range = c(1,6))+
  scale_alpha(range = c(0.3, 0.9))+
  labs(title="Recovered from COVID-19 during the last week",subtitle = "(Top 10 countries)", size= "Recovered cases", alpha = "Recovered cases", colour = "Recovered cases")

ggplot()+
  theme_minimal(base_size=8.5)+
  theme(plot.title = element_text(size = 13,face="bold"),
                plot.subtitle = element_text(size=10))+
  geom_segment(data = recovered_sf_latest_3_change_top10, aes(x= 0,
                                                              xend = change,
                                                              y = reorder(place, change),
                                                              yend = reorder(place, change)),
               colour = "black")+
  geom_point(data = recovered_sf_latest_3_change_top10, aes(x = change, y = reorder(place, change)),
             shape = 21, colour = "white", fill = "firebrick", size= 3)+
  labs(title="Recovered from COVID-19 during the last week",subtitle = "(Top 10 countries)",
       x = "Recovered cases",
       y = "Country")

From graphs we can see that number of recovered cases vary significantly (from 0 to 40000). Europe and North America have most recovered cases, but it would be more meaningful if we examine recovered cases per capita or per cases. The fact that, for example, US and Turkey have the most recovered cases does not means that these countries outperform other countries.

3.5 Population density map for Georgia using classes

ge_adm_0 <- st_read("geo_admbnda_adm0_geostat_20191018.shp")
## Reading layer `geo_admbnda_adm0_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm0_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 13 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 40.00834 ymin: 41.05455 xmax: 46.72672 ymax: 43.5865
## geographic CRS: WGS 84
ge_adm_1 <- st_read("geo_admbnda_adm1_geostat_20191018.shp")
## Reading layer `geo_admbnda_adm1_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm1_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 13 features and 16 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 40.00834 ymin: 41.05455 xmax: 46.72672 ymax: 43.5865
## geographic CRS: WGS 84
ge_adm_2 <- st_read("geo_admbnda_adm2_geostat_20191018.shp")
## Reading layer `geo_admbnda_adm2_geostat_20191018' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\geo_admbnda_adm2_geostat_20191018.shp' using driver `ESRI Shapefile'
## Simple feature collection with 70 features and 19 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 41.54717 ymin: 41.05455 xmax: 46.72672 ymax: 43.25688
## geographic CRS: WGS 84
ge_pop<-read_xlsx("geo_admpop.xlsx",sheet=3)
## New names:
## * `` -> ...8
ge_pop_2 <- ge_pop%>%
  select(ADM2_PCODE, POP_Total)

ge_adm_2_ii <- ge_adm_2 %>%
  select(ADM2_PCODE, ADM2_EN)

ge_adm_2_pop <- left_join(ge_adm_2_ii, ge_pop_2, by="ADM2_PCODE")
ge_adm_1_TBI <- ge_adm_1 %>%
  select(ADM1_PCODE, ADM1_EN) %>%
  rename(ADM_PCODE = ADM1_PCODE,
         ADM_EN = ADM1_EN) %>%
  filter(ADM_EN == "Tbilisi")
ge_adm_1_TBI <- ge_adm_1_TBI %>%
  mutate(POP_Total=1171100)
ge_adm_2_pop <- ge_adm_2_pop %>%
  rename(ADM_PCODE = ADM2_PCODE,
         ADM_EN = ADM2_EN)
ge_adm_2_pop <- rbind(ge_adm_2_pop, ge_adm_1_TBI)
ge_adm_2_pop$area <- st_area(ge_adm_2_pop)
ge_adm_2_pop$area <- set_units(ge_adm_2_pop$area, km^2)
ge_adm_2_pop <- ge_adm_2_pop %>%
  mutate(pop_dens =as.numeric( POP_Total / area))
classes <- classIntervals(ge_adm_2_pop$pop_dens, n = 10, style = "quantile", dig.lab=4)
ge_adm_2_pop <- ge_adm_2_pop %>%
  mutate(percent_class = cut(pop_dens, classes$brks, include.lowest = T, dig.lab = 4))
ggplot()+
  theme_minimal(base_size=8.5)+
  theme(axis.title = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 13,face="bold"))+
  geom_sf(data = ge_adm_0, fill = "grey40", colour = "grey60", size = 0.2)+
  geom_sf(data = ge_adm_2_pop, aes(fill = percent_class), colour = "grey60", size= 0.2)+
scale_fill_brewer(type=seq,palette="RdBu",direction=-1) +
  labs(title = "Population Density in Georgia",
       fill  = "Population/Km2",
        caption = "data: The Humanitarian Data Exchange \nvisual: D.Kajaia")+
  annotation_scale(location = "bl", style = "ticks") +
  annotation_north_arrow(location = "tr", width = unit(0.5, "cm"))

I think population density distribution in Georgia depends on some life conditions. For example, we can see that density is the lowest in the northern regions of Georgia where landscape is hilly and density is the highest in big cities where GDP per capita is very high relative to other regions.

4 Part 4

4.1 Libraries


library(rnaturalearth)
library(tidyverse)
library(tmaptools)
library(lubridate)
library(ggspatial)
library(ggthemes)
library(maptools)
library(GISTools)
library(ggplot2)
library(raster)
library(rgeos)
library(tmap)
library(maps)
library(sf)

4.2 Importing data

setwd("C:/Users/User/Documents/OneDrive/Desktop/R/spatial")
prec <- raster::getData('worldclim', var ='bio', res = 0.5, lon = 42, lat = 44)
gadm_0 <- raster::getData('GADM', country = 'GEO', level = 0)
gadm_0_sf <- st_as_sf(gadm_0)
prec2 <- mask(prec$bio12_17, gadm_0)
roads  <- ne_download(scale = "large", type = "roads", category = "cultural")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\User\AppData\Local\Temp\Rtmp0YCM6x", layer: "ne_10m_roads"
## with 56601 features
## It has 31 fields
## Integer64 fields read as strings:  scalerank question
roads_ge <- gIntersection(roads, gadm_0, byid = TRUE, drop_lower_td = TRUE)
download.file("http://aasa.ut.ee/tbilisi/osm_georgia_selected.zip", destfile="osm_georgia_selected.zip")
unzip("osm_georgia_selected.zip")
file.remove("osm_georgia_selected.zip")
rivers <-  st_read("gis_osm_waterways_free_1.shp")
## Reading layer `gis_osm_waterways_free_1' from data source `C:\Users\User\Documents\OneDrive\Desktop\R\spatial\gis_osm_waterways_free_1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 17123 features and 5 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: 39.52382 ymin: 40.93108 xmax: 46.8273 ymax: 45.74497
## geographic CRS: WGS 84
rivers_2 <- rivers %>% filter(fclass == "river")
rivers_2 <- st_intersection(rivers_2, gadm_0_sf)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
populated_places <- ne_download(scale = "large", type = "populated_places", category = "cultural")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\User\AppData\Local\Temp\Rtmp0YCM6x", layer: "ne_10m_populated_places"
## with 7343 features
## It has 119 fields
## Integer64 fields read as strings:  wof_id ne_id
populated_places_sf <- st_as_sf(populated_places)
populated_places_sf <- populated_places_sf %>% filter(ISO_A2 =="GE")

4.3 Annual Precipitation in Georgia

plot(gadm_0)
plot(prec2, add= T, col = rainbow(20))
plot(rivers_2, add=T,  col="lightblue", lwd = 0.25)
plot(roads_ge, add = T, col = "red", lwd = 0.1)
plot(populated_places_sf, col = "black", add = T)
pointLabel(st_coordinates(populated_places_sf),labels = populated_places_sf$NAME,cex=1)
north.arrow(xb=45.75, yb=43.2, len=0.05, lab="N")
maps::map.scale(x=41.5, y=41.15, relwidth = 0.15,ratio=FALSE, cex=0.9)
title(main ="Annual Precipitation in Georgia", font.main= 1)

4.4 Annual Precipitation in Georgia with ggplot2

re<-as.data.frame(prec2,xy=TRUE,na.rm=TRUE)
str(re)
## 'data.frame':    109577 obs. of  3 variables:
##  $ x       : num  40.2 40.2 40.2 40.3 40.3 ...
##  $ y       : num  43.6 43.6 43.6 43.6 43.6 ...
##  $ bio12_17: int  1265 1223 1200 1215 1226 1223 1279 1319 1259 1285 ...
re <- st_as_sf(re, coords = c("x", "y"), crs = 4326)
roads2<-st_as_sf(roads_ge)
ggplot()+
  theme_grey(base_size=8.5)+
  theme(axis.title = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 13,face="bold"))+
  geom_sf(data = re, aes(color=bio12_17))+
  geom_sf(data=roads2,color="red",size=0.1)+
  geom_sf(data=rivers_2,color="lightblue")+
  geom_sf_label(populated_places_sf,mapping=aes(label=NAME),alpha = 0.6,
                size= 3,expand=TRUE,label.size=0.1)+
  scale_color_gradientn(colours = rainbow(16))+
  labs(color = "mm",
       title = "Annual Precipitation in Georgia",
       caption = "Author: D. Kajaia")+
  annotation_scale(location = "bl", style = "ticks") +
  annotation_north_arrow(location = "tr", width = unit(0.5, "cm"))
## Warning: Ignoring unknown parameters: expand
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

4.5 Average annual precipitation in Georgian regions

gadm_1_sp <- raster::getData('GADM', country = 'GEO', level = 1)
region_Ave_prec <- raster::extract(prec2, gadm_1_sp, fun = mean, na.rm=TRUE, sp = T)
region_Ave_prec_sf <- st_as_sf(region_Ave_prec)
head(region_Ave_prec_sf)
## Simple feature collection with 6 features and 11 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 40.01111 ymin: 41.03851 xmax: 46.72136 ymax: 43.58454
## CRS:            +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
##   GID_0  NAME_0   GID_1       NAME_1 VARNAME_1 NL_NAME_1                 TYPE_1
## 1   GEO Georgia GEO.1_1     Abkhazia   Sokhumi      <NA> Avtonomiuri Respublika
## 5   GEO Georgia GEO.2_1       Ajaria    Batumi      <NA> Avtonomiuri Respublika
## 6   GEO Georgia GEO.3_1        Guria  Ozurgeti      <NA>                 Region
## 7   GEO Georgia GEO.4_1      Imereti   Kutaisi      <NA>                 Region
## 8   GEO Georgia GEO.5_1      Kakheti    Telavi      <NA>                 Region
## 9   GEO Georgia GEO.6_1 Kvemo Kartli   Rustavi      <NA>                 Region
##             ENGTYPE_1 CC_1 HASC_1  bio12_17                       geometry
## 1 Autonomous Republic <NA>  GE.AB 1274.5940 MULTIPOLYGON (((41.9095 42....
## 5 Autonomous Republic <NA>  GE.AJ 1267.8889 MULTIPOLYGON (((41.93469 41...
## 6              Region <NA>  GE.GU 1589.8551 MULTIPOLYGON (((41.72903 42...
## 7              Region <NA>  GE.IM 1089.0719 MULTIPOLYGON (((42.67832 41...
## 8              Region <NA>  GE.KA  708.0507 MULTIPOLYGON (((45.7786 41....
## 9              Region <NA>  GE.KK  623.4208 MULTIPOLYGON (((44.56081 41...
ggplot()+
geom_sf(data = region_Ave_prec_sf, aes(fill = bio12_17), col = "grey", size= 0.25)+
geom_sf_label(data = region_Ave_prec_sf, aes(label= round(bio12_17, 0)), alpha = 0.5)+
scale_fill_gradientn(colours = c("#B0B579","#383b23"))+
labs(title = "Average annual precipitation in Georgian regions", fill = "mm")
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

5 Part 5

5.1 Libraries


library(data.table)
library(tidyverse)
library(lubridate)
library(ggspatial)
library(ggthemes)
library(scico)
library(ggmap)
library(knitr)
library(dplyr)
library(sf)

5.2 Importing GPS data

setwd("C:/Users/User/Documents/OneDrive/Desktop/R/spatial")
url <-  "http://aasa.ut.ee/Rspatial/data/usa_GPS.zip"
download.file(url, "usa_GPS.zip")
unzip("usa_GPS.zip")
file.remove("usa_GPS.zip")
gpsData <- read.csv2("gps_us.csv")

5.3 Cleaning data

## Rows: 185,808
## Columns: 17
## $ user_id        <int> 241, 241, 241, 241, 241, 241, 241, 241, 241, 241, 24...
## $ id             <int> 646380266, 646380267, 646380268, 646380269, 64638027...
## $ restart        <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2...
## $ counter        <dbl> 66718, 66720, 66723, 66730, 66732, 66734, 66736, 667...
## $ time_system    <dbl> 1.524345e+12, 1.524345e+12, 1.524345e+12, 1.524345e+...
## $ time_gps       <dbl> 1.524345e+12, 1.524345e+12, 1.524345e+12, 1.524345e+...
## $ time_system_ts <chr> "2018-04-22 00:10:39", "2018-04-22 00:10:40", "2018-...
## $ time_gps_ts    <chr> "2018-04-22 00:10:36", "2018-04-22 00:10:37", "2018-...
## $ accuracy       <dbl> 12, 16, 16, 6, 4, 6, 6, 6, 4, 4, 6, 6, 6, 6, 6, 6, 6...
## $ altitude       <dbl> -23.991000, -24.298600, -12.608500, -9.616940, -8.17...
## $ bearing        <dbl> 165.957, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ speed          <dbl> 0.104187, NA, 0.000000, 0.000000, 0.000000, 0.000000...
## $ header_id      <int> 226334354, 226334354, 226334354, 226334360, 22633436...
## $ series_id      <int> 1042, 1042, 1042, 1042, 1042, 1042, 1042, 1042, 1042...
## $ year           <int> 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018...
## $ X              <dbl> -121.3842, -121.3842, -121.3842, -121.3842, -121.384...
## $ Y              <dbl> 38.64978, 38.64976, 38.64976, 38.64976, 38.64971, 38...
##        time_system_ts accuracy  altitude bearing    speed         X        Y
## 1 2018-04-22 00:10:39       12 -23.99100 165.957 0.104187 -121.3842 38.64978
## 2 2018-04-22 00:10:40       16 -24.29860      NA       NA -121.3842 38.64976
## 3 2018-04-22 00:10:41       16 -12.60850      NA 0.000000 -121.3842 38.64976
## 4 2018-04-22 00:10:42        6  -9.61694      NA 0.000000 -121.3842 38.64976
## 5 2018-04-22 00:10:43        4  -8.17016      NA 0.000000 -121.3843 38.64971
## 6 2018-04-22 00:10:44        6  -5.67983      NA 0.000000 -121.3843 38.64971

5.4 Importing base map

box <- c(left = min(gpsData2$X)-0.2, bottom = min(gpsData2$Y)-0.2, right = max(gpsData2$X)+0.2, top = max(gpsData2$Y)+0.2)

Ca_map <- get_stamenmap(box, maptype = "terrain", zoom = 9)

5.5 Plot of movement path

ggmap(Ca_map)+
  theme_classic()+
  theme(axis.title = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 13,face="bold"))+
  geom_point(data = gpsData2, aes(y = Y, x=X), colour ="red", size=0.6, alpha=0.5)+
  labs(title ="Movement path on 2018-04-27")+
  annotation_north_arrow(location = "tr", width = unit(0.5, "cm"))

6 Part 6

6.1 Libraries


library(readxl)
library(dplyr)
library(knitr)
library(curl)
library(tmap)
library(sf)

6.2 Importing data

link<-"http://aasa.ut.ee/tbilisi/data/tripadvisor.xlsx"
curl_download(link,destfile="data.xlsx")
data<-read_xlsx("data.xlsx")

6.3 Tidying data

glimpse(data)
## Rows: 987
## Columns: 4
## $ Title   <chr> "Shef's Grandma", "1001 Nights Restaurant", "11 Katkha", "1...
## $ lon     <chr> "44.795958900000002", "44.802948100000002", "44.76982009999...
## $ lat     <chr> "41.695589300000002", "41.693018000000002", "41.72188049999...
## $ reviews <dbl> 35, 8, 1, 97, 45, 82, 28, 2, 2, 16, 2, 4, 194, 81, 81, 1, 2...
data$lon<-as.numeric(data$lon)
## Warning: NAs introduced by coercion
data$lat<-as.numeric(data$lat)
## Warning: NAs introduced by coercion
data2<-data%>%filter(!is.na(lat))%>%filter(!is.na(lon))
dim(data)==dim(data2)
## [1] FALSE  TRUE

6.4 Converting to sf

data_sf<-st_as_sf(data2, coords = c("lon", "lat"), crs = 4326)
glimpse(data_sf)
## Rows: 966
## Columns: 3
## $ Title    <chr> "Shef's Grandma", "1001 Nights Restaurant", "11 Katkha", "...
## $ reviews  <dbl> 35, 8, 1, 97, 45, 82, 28, 2, 2, 16, 2, 4, 194, 81, 81, 1, ...
## $ geometry <POINT [°]> POINT (44.79596 41.69559), POINT (44.80295 41.69302)...

6.5 Map of restaurants in Tbilisi

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(data_sf)+
  tm_dots(size = "reviews", col = "red", alpha = 0.5, scale = 1, border.col="black", sizes.legend=c(0.5,2))
## Legend for symbol sizes not available in view mode.

7 Part 7

7.1 Libraries


library(rnaturalearth)
library(tidyverse)
library(wbstats)
library(dplyr)
library(tmap)

7.2 Choosing indicator

wbstats_indicators <- wb_cachelist$indicators
write.csv(wbstats_indicators, file = "wbstats_indicators.csv")

After exporting indicators to CSV files I chose “Renewable electricity (% in total electricity output): Share of electrity generated by renewable power plants in total electricity” with ID - “4.1_SHARE.RE.IN.ELECTRICITY”

data<- wb(indicator ="4.1_SHARE.RE.IN.ELECTRICITY",
                   mrv = 1, return_wide = TRUE)
glimpse(data)
## Rows: 232
## Columns: 5
## $ iso3c                         <chr> "ABW", "AFG", "AGO", "ALB", "AND", "A...
## $ date                          <chr> "2015", "2015", "2015", "2015", "2015...
## $ iso2c                         <chr> "AW", "AF", "AO", "AL", "AD", "AE", "...
## $ country                       <chr> "Aruba", "Afghanistan", "Angola", "Al...
## $ `4.1_SHARE.RE.IN.ELECTRICITY` <dbl> 14.8561614, 86.0501113, 53.1749283, 1...
colnames(data)[5]<-"Value"
glimpse(data)
## Rows: 232
## Columns: 5
## $ iso3c   <chr> "ABW", "AFG", "AGO", "ALB", "AND", "ARE", "ARG", "ARM", "AS...
## $ date    <chr> "2015", "2015", "2015", "2015", "2015", "2015", "2015", "20...
## $ iso2c   <chr> "AW", "AF", "AO", "AL", "AD", "AE", "AR", "AM", "AS", "AG",...
## $ country <chr> "Aruba", "Afghanistan", "Angola", "Albania", "Andorra", "Un...
## $ Value   <dbl> 14.8561614, 86.0501113, 53.1749283, 100.0000000, 86.1167002...

7.3 Downloading map

countries <- ne_download(scale = 50, type = 'countries', category = 'cultural', returnclass = "sf")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\User\AppData\Local\Temp\Rtmp0YCM6x", layer: "ne_50m_admin_0_countries"
## with 241 features
## It has 94 fields
## Integer64 fields read as strings:  POP_EST NE_ID
countries2 <- countries %>%
  dplyr::select(ADM0_A3, NAME, GDP_MD_EST, POP_EST)

7.4 Joining data table with geospatial layer

data_joined <- left_join(countries2, data, by = c("ADM0_A3" = "iso3c"))

7.5 Creating plots

tmap_mode("plot")
## tmap mode set to plotting
tm_shape(data_joined)+
  tm_polygons(col = "Value",
              style = "kmeans",
              n=5,
              palette=c('#d7191c','#fdae61','#ffffbf','#a6d96a','#1a9641'),
              border.col = "grey50",
              lwd = 0.2,
              title = "Percentage",
              colorNA  = "grey80")+
  tm_style("natural")+
  tm_format("NLD")+
  tm_layout(main.title = "The Share of Renewable Electricity",
            legend.bg.color = "white",
            legend.bg.alpha = 0.6)+
  tm_credits("Data: World Bank\nMap: D. Kajaia",
             bg.color = "white",
             bg.alpha = 0.7)

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(data_joined)+
  tm_polygons(col = "Value",
              style = "kmeans",
              n=5,
              palette=c('#d7191c','#fdae61','#ffffbf','#a6d96a','#1a9641'),
              border.col = "grey50",
              lwd = 0.2,
              title = "Percentage",
              colorNA  = "grey80")+
  tm_style("natural")+
  tm_format("NLD")+
  tm_layout(main.title = "The Share of Renewable Electricity",
            legend.bg.color = "white",
            legend.bg.alpha = 0.6)+
  tm_credits("Data: World Bank\nMap: D. Kajaia",
             bg.color = "white",
             bg.alpha = 0.7)
## Credits not supported in view mode.
## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs).
## Need '+proj=longlat +datum=WGS84'